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ABSTRACT 

In this paper, globular star clusters which contain a sub-system of stellar-mass 
black holes (BH) are investigated. This is done by considering two-component models, 
as these are the simplest approximation of more realistic multi-mass systems, where 
one component represents the BH population and the other represents all the other 
stars. These systems are found to undergo a long phase of evolution where the centre of 
the system is dominated by a dense BH sub-system. After mass segregation has driven 
most of the BH into a compact sub-system, the evolution of the BH sub-system is found 
to be influenced by the cluster in which it is contained. The BH sub-system evolves in 
such a way as to satisfy the energy demands of the whole cluster, just as the core of a 
one component system must satisfy the energy demands of the whole cluster. The BH 
sub-system is found to exist for a significant amount of time. It takes approximately 
lOtrh.i, where iih,i is the initial half-mass relaxation time, from the formation of the 
compact BH sub-system up until the time when 90% of the sub-system total mass 
is lost (which is of order 10^ times the half-mass relaxation time of the BH sub- 
system at its time of formation). Based on theoretical arguments the rate of mass 
loss from the BH sub-system (M2) is predicted to be — /3CM/(airh), where M is the 
total mass, i^h is the half-mass relaxation time, and a, /3, ^ are three dimcnsionless 
parameters (see Section [2] for details) . An interesting consequence of this is that the 
rate of mass loss from the BH sub-system is approximately independent of the stellar 
mass ratio {m2/mi) and the total mass ratio (M2/M1) (in the range m2/mi > 10 and 
M2/M1 ~ 10~^, where mi, m2 are the masses of individual low-mass and high-mass 
particles respectively, and Mi, M2 are the corresponding total masses). The theory is 
found to be in reasonable agreement with most of the results of a series of N-body 
simulations, and with all of the models if the value of C is suitable adjusted. Predictions 
based on theoretical arguments are also made about the structure of BH sub-systems. 
Other aspects of the evolution are also considered such as the conditions for the onset 
of gravothermal oscillation. 

Key words: globular clusters: general; methods: numerical; methods: n-body simu- 
lations. 



1 INTRODUCTION 

Hundreds of stellar mass black holes (BH) can form within 
a massive globular cluster (see Kulkarni, Hut & McMillan 



( 1993 1, Sigurdsson & Hernquist ( 1993 1 and Portegies Zwart 



fc McMillan (20001). Some of the BH might escape at the 



time of their formation due to large natal kicks. However the 



subject of natal kicks for BH is still under debate ( Repetto 



2012 1 and it is possible that the largest BH may form with- 
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out any supernova explosion (Fryer 19991. Uncertainty in 



the natal kicks leads to uncertainty in the initial size of the 
BH population. As the BH are more massive than the other 
stars in the system, any retained BH will undergo mass seg- 
regation and almost all are likely to become concentrated 
in the centre of the system, eventually forming a compact 
sub-system. 

The mass of the BH sub-system decreases over time 
because BH binaries form in the dense core of the BH sub- 
system, causing the ejection of single BH and ultimately the 
binaries themselves through super-elastic encounters (see 



Kulkarni, Hut & McMillan (19931, Sigurdsson & Hernquist 
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(19931 



(20101 



Portegies Zwart & McMillan' ("2000 1 



Banerjee et al 



Downing et al| (|2010[), [Aarsetli (2012l). Early work 



by Kulkarni, Hut fc McMillan (19931 and 



Sigurdsson & 



Hernquist ( 1993 1 seemed to indicate that the BH popula- 



tion will become depleted over a relatively short timescale. 
This conclusion was reached in part by treating the BH sub- 
system as if it were an independent system once most of the 
BH had segregated to the centre of the cluster. 



Merritt et al (20041 and Mackey et al (20081 found that 



heating by a retained population of BH causes large-scale 
core expansion in massive star clusters. They suggest this 
may partly explain the core radius-age trend observed for 
such objects in the Magellanic Clouds. The BH binaries that 
are formed in the core of the BH sub-system are an inter- 
esting class of objects in their own right, especially as the 
merger of two BH may be detectable as a source of grav- 



itational waves ( Portegies Zwart & McMillan ( 2000 1 and 



Banerjee et al (20101). It has even been suggested that star 
clusters consisting almost entirely of BH, known as dark star 



clusters, could exist (Banerjee & Kroupa 20111. Dark star 



clusters could be created if the stars in the outer parts of 
a larger system were stripped away by a strong tidal field, 
leaving behind the BH sub-system. If one were to observe 
the few remaining stars in these systems they would appear 
to be super virial, as the velocity dispersion of the remaining 
stars would be enhanced by the unseen BH. 

[Breen fc Heggie] ( |2012a[ ) investigated the evolution of 
two-component models and found that, within the parame- 
ter space they considered, the stability of the two-component 
system against gravothermal oscillations was dominated by 
the heavy component. They only considered systems with a 
total mass ratio of order M2/M1 > 10"\ where M2 (Mi) is 
the total mass of the heavy (light) component. However the 
mass ratio of a system containing a BH s ub-system would 
only be expected to be M2/M1 ~ 10~^ (Portegies Zwart 



|fc McMilla"n||2000| l, where M2 is the total mass of the BH 
sub-system, which is smaller by an order of magnitude than 
any of the systems studied by Breen & Heggie (2012a). As 



Henon's Principle ( |Henon|1975[ ) states that the energy gen- 
erating rate of the core is regulated by the bulk of the sys- 
tem, it seems unlikely that the approach of [Breen fc HeggTe] 



( [2012a| ) is appropriate in this case due to the small value of 
M2/M1. 

This paper is structured as follows. In Section |2] some 
theoretical results are derived and discussed. This is followed 
by Section|3J where the theoretical results regarding the size 
of the BH sub-system are tested using both gas models and 
direct N-body runs. Section |4| contains empirical results re- 
garding the mass loss rates from BH sub-systems and a 
comparison between the empirical results and the theory 
of Section [2] The qualitative behaviour of these systems is 
also discussed in this section. Section [S] is concerned with 
gravothermal oscillations in systems containing a BH sub- 
system. Finally Section |6| consists of the conclusions and a 
discussion. 



2 THEORETICAL UNDERSTANDING 

2.1 BH sub-system: half- mass radius 

Here we will consider aspects of the dynamics of a system 
containing a BH sub-system. We will assume that the system 



is Spitzer unstable (ISpi tzer|1987 l and that the total mass of 
the BH sub-system (M2) is very small compared to the total 
mass of the system (M). (Since the Spitzer stability criterion 
is (M2/Mi)(m2/mi) 2 < 0.16, these assumptions are consis- 
tent provided the stellar mass ratio -012/ mi is large enough). 
We will also assume that the initial state of the system has 
a constant mass density ratio between the two components 
throughout and that the velocity dispersions of both com- 
ponents are equal at all locations. If this is the case then the 
system would first experience a mass segregation-dominated 
phase of evolution which lasts of order {mi/m2)tcc (Fregeau 



[et al|20(J2[ and references therein), where 1712 (rni) is the stel- 
lar mass of the BH (other stars), and ice is the core collapse 
time in a single component system, although technically for 
the outermost BH mass segregation can last much longer 
than {mi/m2)tcc (see Appendix [b] for details). 

If we consider the 50% Lagrangian shell of the heavy 
component, initially it will be approximately the same size 
as the 50% Lagrangian shell of the entire system. As the 
BH lose energy to the other stars in the system the 50% 
Lagrangian shell of the BH component contracts. The shell 
will continue to contract until the energy loss to the light 
component is balanced by the energy the shell receives from 
the inner parts of the BH component. As we have assumed 
that the system is Spitzer unstable, it follows that a temper- 
ature difference must remain between the two components 
and thus there is still a transfer of energy between the two 
components. As the total mass of BH is small, the contrac- 
tion of the 50% Lagrangian shell of the heavy component 
continues until the system is concentrated in a small region 
in the centre of the system. This is what we call the BH sub- 
system. The BH sub-system is very compact and therefore 
rapidly undergoes core collapse. The subsequent generation 
of energy by the formation of BH binaries, and interactions 
of BH binaries with single BH, support the 50% Lagrangian 
shell of the heavy component. 

In the present paper, we will assume that the main 
pathway for the transport of thermal energy throughout the 
system is as follows: energy is generated in the core of the 
BH sub-system (we are assuming that the BH core radius 
is much smaller than rh,2, the half mass radius of the BH 
sub-system), then the energy is conducted throughout the 
BH sub-system via two-body relaxation just as in the con- 
ventional picture of post-collapse evolution; but in the stan- 
dard one-component setting this flux causes expansion and, 
ultimately, dispersal of the system. In a two-component sys- 
tem, however, the coupling to the lighter component changes 
the picture dramatically. We will assume that at a radius 
comparable with rh,2 most of the energy flux is transferred 
into the light component, where it then spreads throughout 
the bulk of the light system. These assumptions will hold if 
most of the heating, either direct (heating by reaction prod- 
ucts which remain in the cluster) or indirect, initially occurs 
within the BH sub-system rather than within the regions 
dominated by the light component (see Appendix [A] for a 
discussion on this issue). A similar assumption is made in 



one-component gas models (Goodman 1987 Heggie fc Ra- 
mamani [1989 1 which have been shown to be in good agree- 
ment with direct N-body models ( [Bettwieser fc Sugimoto' 
19851). 



From Henon's Principle (Henon 1975 1 we argue that the 



rate of energy generation is regulated by the energy demands 
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of the bulk of the system. For the systems we are considering 
here the bulk of the system is in the hght component as 
Ml ^ 0.99M for M2/M1 = 10"^. The energy demands of a 
system are normaUy thought of as the energy flux at the half- 
mass radius, which is of order |i5|/frh, where E is the total 
energy of the system and trh is the half mass relaxation time. 
Under our assumptions the energy flux must be supplied by 
the BH sub-system, which ultimately must be generated by 
binaries in the core of the BH sub-system. However, we can 
ignore the details of how the energy is actually generated and 
consider the BH sub-system itself as the energy source for 
the cluster as a whole. In this picture, the energy exchange 
(_Eex) between the BH sub-system and the light component 
must balance the flux at the half mass radius of the light 
component i.e. 






\E,: 



The BH sub-system is concentrated in the centre of the 
cluster, and therefore the half-mass relaxation time in the 
BH sub-system (trh,2) is quite short. It follows that the flux 
at the half-mass radius of the BH sub-system is quite high. 
The flux at the half-mass radius of the BH sub-system (rh,2) 
is of order \E2\/tih,2, where E2 is the energy of the BH sub- 
system. Most of the energy that passes rh,2 must be trans- 
ferred to the light component or else the BH sub-system 
would rapidly expand until the flux around rh,2 is compa- 
rable to the rate of energy exchange. Therefore the energy 
exchange rate between the two components must be approx- 
imately equal to the flux at the half mass radius of the heavy 
sub-system i.e. 



irh,2 



\E,: 



All this leads to the conclusion that the flux at r^ must 
be balanced by the flux at rh,2 i.e. 



which can be rearranged as 

\E\ 



irh,2 



irh,2 



(1) 



Using the definition of t^h as given in |Spitzer] ( |1987[ ) (and an 
equivalent definition for frh,2), the right hand side of equa- 
tion [l] becomes 



(2) 



where m is the mean mass, and In A, lnA2 are the coulomb 
logarithms of the entirely system and the BH sub-system 
respectively. 

Using equation [l] the left hand side of equation [2] be- 
comes 



irh 


^ Nh^m^ InAa 


M 2 r^ m,2 In A2 


irh,2 


-13 1 

N^r^^mi In A 


1 3 



\E\ 
\Eh\ 



Ma^ 



12^2 



MVh,; 



(3) 



where we have estimated the squared one dimensional ve- 
locity dispersions, 0-2 and a^, by assuming that both the 
system and the sub-system are in virial equilibrium, so that 



o-^ ~ 0.2GM/rh and a^ ~ 0.2GM2/rh,2. Putting the above 
equations together we have 

1 s 
M^rh,2 M2m2r^lnA2 

and then by rearranging this we have 



M^ m2 In A2 
]\/[l m InA 



(4) 



This result implies that for a fixed total mass ra- 
tio M2/M1 (« M2/M) and ignoring the variation of the 
coulomb logarithms, the ratio of rh,2/»'h grows with increas- 
ing stellar mass ratio m2/m.. 



2.2 BH sub-system: core radius 



In Section 2.1 rh,2/rh was estimated by assuming that the 



energy flow in the BH sub-system balances the energy flow 
in the bulk of the system (i.e. the other stars). In order for 
equation HI to hold it is assumed that the BH core radius 
rc,2 must be significantly smaller than rh,2 and that the BH 
sub-system is actually capable of producing and supplying 
the energy required by the system. Usually it is assumed 
that the core (in this case the core of the BH sub-system) 
adjusts to provide the energy required. In this section we 
estimate the size of the core, and use the estimate to place 
a condition on the validity of our assumptions. 

As most of the mass within rh,2 is BH, we may treat 
the BH sub-system as a one-component system and we can 
make use of the standard treatments of one-component sys- 
tems. In balanced evolution (i.e. a situation in which en- 
ergy is produced at the rate at which it fiows over the 
half-mass radius), the rate of energy production is given 
by _E = (|-E2|/irh,2)C2, where C,2 is a constant (for a one- 
component model ("2 ~ 0.0926, see Henon (19651). We will 



follow the derivation in Heggie & Hut (20031 page 265) of 



the dependence of rc/rh on A'^ for a one component model, 
although here it will be necessary to keep track of the numer- 
ical constants and account for the fact that the properties 
correspond to those of the BH sub-system. In order to de- 
rive a condition on re, 2/^11,2 it is necessary to express E, E2 
and trh,2 in terms of the other properties of the BH sub- 
system. E « Mc,2e where Mc,2 is the BH core mass and e 
is the energy generating rate per unit mass (Heggie fc Hut] 
|2003 |). Mc,2 and e can be expressed in terms of Pc,2 (the 
central mass density of BH), re, 2 and crc,2 (the central one- 
dimensional velocity dispersion of the BH), which results in 



E = 85 

'^c,2 

As in the previous section we shall use 

ii2 ~ U./ . 

rh,2 

It will be convenient to use a different but equivalent defini- 
tion of trh,2 ( |Spitzer||1987| ) rather than the one used in the 
previous section, i.e. 



irh,2 = 



0.195cr^ 



Gm2Ph,2lnA2 
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where (T2 is the one-dimensional velocity dispersion of the 
BH inside rh,2 and ph,2 is the mean mass density of the BH 
inside rh,2- 

Using E — (|-B2l/trh.2)C2 all of the above equations can 
be combined into 

83G m2Pc,2'^2rc,2rh,2 ~ M2 cTc, 2Ph,2C2 In A2. 



This can be simplified by using inGpca^c 



9CT(, 2- Also the 



BH sub-system is expected to be nearly isothermal inside 
rh,2', therefore 0-2 ~ crc,2 and furthermore it follows that 
P2 oc r~^ inside rh,2- Taking all of this into account and by 
rearranging the above we have 



re, 2 

fh,2 



A^o" 



VC2lnA2/ 



(5) 



2.3 Limitations of the theory 

One of our assumptions was that re, 2 -C rh,2 and now we 
can derive an approximate condition for the validity of the 
theory. As A'^2 is small we shall take lnA2 ~ 1. As the 
entire system is in balanced evolution it is also true that 
{\E2\/trh,2)C2 = (|-E|/trh)C! where C is a dimensionless pa- 
rameter defined implicitly by the equation E — (|_E|/£rh)Cj 
we expect |-E2|/irh,2 ~ |^|/irh (equation [l| and so for the 
purposes of our estimate we can assume that (2 ~ C- There- 
fore rc,2 ^ rh,2 provided A^2 ^ 40. This value is only a rough 
guide, and what is important to take from this result is that 
for sufficiently small N2 the theory in Section |2.1| breaks 
down. 

As M2 decreases the BH sub-system will ultimately 
reach a point where it can no longer solely power the expan- 
sion of the system by the mechanism we have considered (i.e. 
formation, hardening and ejection of BH binaries by inter- 
action amongst the BH sub-system). After this point it may 
be possible for BH binaries to generate the required energy 
through strong interaction with the light stars. However this 
will probably require a much higher central mass density of 
the light component than at the time of formation of the 
BH sub-system, as at this central mass density interactions 
between the light stars and the BH binaries are expected to 
be much less efficient at generating energy than interactions 
between single BH and BH binaries. This implies a poten- 
tially significant adjustment phase towards the end of the 
life of the BH sub-system, as is illustrated by an N-body 
model in Section HI (see Fig. ml . 

In Section |2.l| we made the assumption that the BH 
sub-system was Spitzer stable. However, as pointed out to 
us by Sambaran Banerjee (private communication), it is 
also possible that as M2 decreases a point may be reached 
were the sub-system becomes Spitzer stable. If so, the two 
components could reach equipartition at the centre, i.e. 
wi2cr2 = miai, and in that case our assumption that heat 
fiows from the heavy component to the light component is 
false. The temperature ratio of the two component may be 
estimated by 



7712(72 


m2 M2 Th 


/ m2 N 1 / M2 \ f / In A2 
'^KmiJ \Mi) V InA 


miaf 


nil Ml rh,2 



were we have made u se o f equation HI and the assump- 
tions made in Section 2.1 Initially (m2(T2)/(micri) > 1 



{{m2(J2) / {niiaf) = 1) as BH escape and M2 decreases. 
By setting {m20-2) / [niiai) — 1, ignoring the variation of 
coulomb logarithms and taking each side to the power 2/5 
we arrive at 



M2 
Ml 



Vmi / 



where C is a constant. This is exactly the same form as the 
Spitzer stability criterion ( Spitzer] 1987 1, the only difference 
being that we have not specified the constant C. Again it 
follows that the theory of Sections |2.1| and |2.2| will fail when 
M2 becomes too small, and that the limiting value of A/2 is 
smaller for larger m2/mi. 

We will now briefly consider the case of a Spitzer stable 
BH sub-system. As the BH move more slowly than the other 
stars they still concentrate in the centre of the system. If the 
heavy component still dominates within rh,2 then the BH 
sub-system is self-gravitating and 



2 
0-2 



GM2 
rh.2 



(6) 



still holds. From equipartition of kinetic energy we have 

mi 2 "T-i GM GM2 

CTl ■■ ~ . 

"12 m2 Th rh,2 

This result can be rearranged as 

rh,2 7712 A/2 

Th Till M 

In equation [6] there is a difi'erent dependence of rh,2/rh on 
7712/1111 and A/2 /A/ than in equation [s] 



2.4 Evaporation rate 

We will now consider the evaporation rate as a result of 
two-body encounters for the BH sub-system. It is important 
to note that evaporation is only one of the mechanisms by 
which BH are removed from the system. Another important 
mechanism as already discussed is ejection via encounters 
involving BH binaries and single BH. In this section we will 
ignore this effect although it will be considered in detail in 
the next section. 

The one dimensional velocity dispersion of the BH sub- 
system has the following dependence on A/2 and rh,2 (as- 
suming it is nearly self gravitating) : 



2 
0-2 



GM2 

rh,2 



From the previous section 



know that 



V rh / 



M^ ni2 In A2 



Therefore, if we consider the post 



Mi m InA 

collapse evolution of a series of models with different values 
of 1712/711, at the same values of rh, A/2 and M, as 1112/711 
increases so does rh,2 and therefore o"! decreases. Here we 
have ignored the variation of the coulomb logarithms; for a 
system with A*' = 10*^ and A/2/A/1 = 10"^, the variation of 
lnA2/lnA is a factor of 2.2 between 1112/7111 = 10 and 50 if 
A2 = 0.02A'2 and A = 0.02N. The source of the value 0.02 
is 



Giersz fc Heggie (19961 



but for fixed 1112/7111 it decreases towards equipartition 



The mean-square escape velocity is related to the mean- 
square velocity of t he system (see Spitzer ( 1987 1 and Binney 
& Tremaine (2008 1 ) by v'^ = 4d'^. The mean-square velocity 
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of the system is dominated by the hght component and re- 
mains approximately fixed with varying m^/m. Therefore, 
as m-ilm increases the mean-square velocity of the BH sub- 
system decreases relative to the mean-square escape veloc- 
ity. This implies that systems with higher m2/m lose a lower 
fraction of their stars by evaporation per irh,2; in fact we will 
show in the next paragraph that escape via evaporation is 
negligible. 

A rough estimate of the fraction of stars lost by evapora- 
tion each trh can be calculated from the Maxwellian velocity 
distribution ( |Spitzer||1987[ [Binney fc Tremain"e "2008). This 
is done by assuming that the fraction of stars with v greater 
than Vc in a Maxwellian velocity distribution is removed each 
trh, i.e. 



dN 
dt 



N 

irh 



7, 



where 7 denotes the fraction of stars removed. Its value for 
a one component model is 7 = 7.38 x 10~^. In order to es- 
timate the value of 7 for a two component model we need 
to know the relationship between «! and v\, where wf (vi) 
is the mean-square velocity of the heavy (light) component. 
If the system was Spitzer stable then m,2W2 = inivf; how- 
ever the systems we are considering are not Spitzer stable 
because of the large stellar mass ratios, and therefore it is 
expected that m2V2 > m\v\. Over a period of time where 
M2 and rh remain roughly constant then V2 and vl will be 
approximately constant. Therefore over the same time pe- 
riod we have 7TI2W2 ~ kmiv\, where fc is a constant. Using 



a two-component gas model (see Heggie & Aarseth ( 1992 1 



and Breen & Heggie (2012al), for the range of parameters 
in Section |3j fc was found to be ^ 2. Assuming a stellar mass 
ratio of 10 and letting fc = 2 to insure the highest possible 
value of 7 leads to 7 = 5.87 x 10"^'^. This exceedingly small 
value of 7 is a result of the fact that the Maxwellian velocity 
distribution drops exponentially with increasing velocity, so 
that even a slight increase in escape velocity leads to a much 
smaller value of 7. 

Based on this approximate theory we can conclude that 
mass loss from evaporation due to two body encounters is 
not significant for the case of BH sub-systems. It is worth 
noting (based on the arguments in this section) that con- 
straints based on evaporation timescales (for example see 



Maoz 1998 1 which only take into consideration the poten- 



tial of the BH sub-system are not generally valid if the sub- 
system is embedded in a much more massive system. BH 
which escape the sub-system in two-body encounters gener- 
ally cannot escape from the deep potential well of the sur- 
rounding system. Instead, they return to the sub-system on 
the mass segregation/dynamical friction timescale. 



2.5 Ejection rate 

Dynamical evolution of BH binaries and ejection of BH is 
an energy source which is assumed in the present paper to 
com ply with Henon's Principle. As has been stated in Sec- 
tion 



2.1 



Henon's Principle states that E is regulated by the 



energy demands of the bulk of the system, i.e. 



E: 



\E\ 



where (^ is a constant. For systems with M2 <^ M, (|_B|/frh)C 
is determined mainly by the properties of the light com- 
ponent and is approximately independent of M2 and m2. 
Therefore the energy generation rate is also approximately 
independent of M2 and m2 . 

The encounters which generate energy (either by forma- 
tion of binaries or their subsequent harding) happen where 
the density is highest, in the core of the BH sub-system. 
As the BH are concentrated in the centre of the system, 
through mass segregation, we may assume that encounters 
which generate energy predominantly occur between BH bi- 
naries and single BH. The BH sub-systems considered in 
this paper consist of BH with identical stellar mass, there- 
fore the mechanism by which energy is generated in the BH 
sub-systems is similar to that for a one-component system. 
The two key differences for the BH sub-system are that the 
escape potential is elevated and the size of the system is reg- 
ulated by the much more massive system of light stars (see 
equation |4| . 

In a one-component system each hard binary formed 
in the core on average contributes a fixed amount of en- 
ergy (X m(j>c (where (pc is the central potential) before be- 
ing ejected from the system (see for example Heggie & Hut 



2003 1 . Typical estimates of the average energy each hard bi- 
nary contributes in a one compone nt system are ~ 7.5m(l>c 
( Goodman] fl984| and si 8.27mcf)c ( [Heggie fc Hut| [2b03p] 



Also on average each hard binary causes the ejection of a 



fixed number of stars. Goodman ( 1984 1 estimated this to be 



approximately 6 stars (including the binary itself) and Heg- 



gie fc Hutj (2003) estimated this to be approximately 3 stars 
(excluding the binary itself). The situation is similar for a 
BH sub-system and we can assume that the mass ejected and 
the average contribution per hard BH binary is the same as 
for the one-component case. Furthermore as mass loss due 
to evaporation is negligible for a BH sub-system (see Section 
2.4 1, the loss of mass from the sub-system is always associ- 
ated with energy generation. Therefore we can express the 
rate of energy generation in the core in terms of mass loss. 



E « /3M2 0C 



(8) 



where /3 is a constant; /? ~ 2.2 in the one component case, 
where we have used the values of energy generated and mass 



lost given in Heggie fc Hut (20031, adjusted to account for 



the energy generated and mass lost in the escape of the bi- 
nary itself, « lO.6m,20c and « 4.7m2 receptively. Since E is 
regulated by the light component (equation [tI, we can use 
equation [S] to estimate the rate of mass loss. Note that the 
estimates in this paragraph are entirely theoretical, with- 
out detailed numerical support especially for the value of /3. 
Note also that the estimate ignores the heating effect of en- 
counters which do not lead to ejection once the binary has 
reached a sufficient binding energy for ejection to be likely. 
We will now show that (f>c is approximately indepen- 
dent of the properties of the BH sub-system. This will be 
done by showing that the main contribution to the central 
potential is from the light component. We can estimate the 
contribution of the lights to ^c to be <j!>i « — GM/rh and 



(7) 



^ Note there is an error in [Heggie &: Hut| ( |2003| l p. 225: the con- 
stant is stated incorrectly but the correct value can be obtained 
by evaluating the formula given on the same page. 



© 2012 RAS, MNRAS 000,[lH20l 



6 P. G. Breen and D. C. Heggie 



the contribution of the BH to be 



regime of interest M2/M = 10 ^ and rh,2/rh = 10 ^ (for 
typical values oirh,2/rh- see Tablell]). Therefore, 



-GA/2/7'h,2- In the and if we integrate the above equation we get 



(h ^ M2 Th 
■^1 ~ M rh,2 



10" 



and we can approximate <j>c by (jii (see also Section 3.3 1. 

We can now use E « /3Af2(/'c to make an estimate of 
the mass loss rate: from equations [7] and |8] we have 



pM2^. ~ ^C, 

Jrh 



and so 



M2 ~M 



\E\ 1 C 



The term E/{M(j)c) is dimensionless and approximately in- 
dependent of the properties of the BH sub-system; we wiU 
use a to represent its value. For a Plummer model a « 0.15, 
however during core collapse |0c| increases while E and M 
remain approximately constant, resulting in smaller values 
of a. The two-component gas models used in Section [S] in- 
dicate a value of a « 0.13. We now have 



M2 



M_aC 



Scaling to the values oi a, ( and /3, we have 



M2 



-0.0061 



M 



C 2.2 



"irh 0.15 0.09 /3 



(9) 



(10) 



By this estimate the sub-system should last ~ 1.6irh to 
3.3frh, for M2/M1 = 0.01 to 0.02 and canonical values of a, 
P and ^. The important point to take from this result is that 
the rate of mass loss from the BH sub-system depends on 
the half-mass relaxation time of the whole system and not 
on any property of the BH sub-system. 

While a system is in balanced evolution, the only pa- 
rameter that varies significantly (over a timescale where M 
is negligible) in the right hand side of the above equation 
is tih, due to the increase in rh (Here we assume that the 
system is isolated; the case of a tidally limited system is 
considered in the following section.). Therefore, for a par- 
ticular system equation [lO] can be expressed in the form 

where C is a constant (C — -r. — : '- where 



M2 



-Crl 



/3 fA,i 

aC,l P ~ 6.1 X 10~^ for canonical values of a, /3, and C, and 
i denotes values at the start of the balanced evolution). r\ 
itself is a function of time which can be derived from the 
relation Vh/rh ~ C/trh, which follows in turn from equa- 
tion I7I if we assume E ex GM^/r-ti and we assume mass 

3 
loss from the entire system is negligible. Since frh oc r^ it 

3C 



follows that Th i; rh 



V 2trh,i 



tc 



c)) , where powered 

expansion starts at time tec (the reason for this notation will 
become clear in Section U|. Therefore M2 can be expressed 
as 



M2 



Crl 



1 + 



3C 

2trh.: 



{t - tec) 



(11) 



M2 ~ M2 



2a 



Mln 



(^ + ^(*-^-))- 



(12) 



Throughout this section it has been assumed that ^ is 
a constant. However, as discussed in Section |4J C, has been 
found to vary with time in situations where the BH sub- 
system cannot provide enough energy for balanced evolu- 
tion. Even if C, is time-dependent equation |9] and the equa- 
tion 



1 . 

Th 



trh 



c, 



are still expected to hold (under the other assumptions made 
in this section). These two equations can be combined into a 
single equation which relates M2 to rh and has no explicit (" 

a 1 
dependence (i.e. M2 = —M— — rh). The resulting equation 

p rh 
can be easily solved (assuming that the variation of M and 
a//3 are neglected) and its solution is 



M2 =M2 



M^ln^ 
P rh,i 



(13) 



This result implies that, regardless of C,, systems with the 
same M2,i and rh,i should evolves along the same curve in 
M2, rh space. 



2.6 Tidally limited systems 

In this section we will briefly consider the theory of tidally 
limited systems containing BH sub-systems. In Henon's 



tidally limited model (Henon 19611, the rate of mass loss 
is 



M 



M 

trh 



c 



where ^ is a constant (^ = 0.045). In Section 2.5 
of the same form was found for M2, i.e. 



(14) 
an equation 



M2 



A£aC 



The relation between M and M2 can be found by simply 
dividing these two equations, which results in 



M2 
M 






0.11- 



C 2.2 0.045 



0.15 0.0725 /3 C 



Therefore M2/M is a constant. For canonical values of 
the constants in the above equation M2/M ~ 0.11. Note 
that the tidally limited model has a different value of ^ 
{C, ~ 0.0725, sec Hcnon '1961) than for an isolated model 
(C = 0.0926, sec Henon 1965 ). The constant value of M2/M 
implies that for two-component systems there is a threshold 
value of M2/M1 at ~ 10"^ above which M2/M (and hence 
M2/M1) is expected to grow with time and below which 
M2/M decreases with time. In other words if A/2 /Mi > 10"^ 
then the system is expected to become more BH dominated. 



ultimately becoming a so-called dark star cluster (Baner 



jee fc Kroupa 20111. Alternatively if M2/M1 < 10"' then 



the BH sub-system is expected to dissolve. In Section |6.3[ 
where two-component parameter space is classified into dif- 
ferent regions (see Fig. |16[ ), there is already a distinction at 
roughly M2/M1 ~ 10~ (between region II & HI systems). 
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Table 1. Values of ?'ii,2/''h in post collapse evolution {N = 32fe). 
These values where measured over li^ji i after a time of at least 
2tcc, where tec is the time of core collapse. For low AI2/M1 {< 0.1) 
there is a clear trend of increasing r^^ 2/''h with increasing m2/m\ . 
The total mass ratios of 0.5 and 0.1 have also been included to 
demonstrate that there is an inverse dependence of r^ 2/''h on 
m2/m\ for the largest value of M2/M1 considered. See text for 
more details. 



3 

>^ 

-a 
o 



m2 \ M2 
mi ^ All 


0.5 


0.1 


0.05 


0.02 


0.01 


-a 


100 


0.37 


0.34 


0.29 


0.24 


0.21 




50 


0.38 


0.27 


0.24 


0.19 


0.18 




20 


0.40 


0.21 


0.18 


0.15 


0.13 




10 


0.42 


0.18 


0.14 


0.11 


0.09 




5 


0.44 


0.18 


0.11 


0.07 


0.04 




2 


0.50 


0.20 


0.10 


0.03 


0.02 






1.5 2 2.5 

Time (t^^^ ;) 



based on other reasons discussed in that section. The the- 
ory in this section can be viewed as another reason for the 
distinction. 

It is important to note that this result has yet to be rig- 
orously tested because tidally hmited systems are not con- 
sidered further in this paper and the exact threshold value 
is likely to depend on a number of astrophysical issues (e.g. 
initial mass function, tidal shocks etc). Indeed while equa- 
tion [14] is a reasonable approximation, |Baunigardt| ( |2001[ ) 
showed that the time scale of escape depends on both irh 
and the crossing time. Nevertheless it seem likely that a 
threshold value of AI2/M1 exists even for more realistic sys- 
tems, although it may have some dependence on the other 
properties of the system (e.g. 7712 /mi). 



3 DEPENDENCE OF rh,2/rh ON CLUSTER 
PARAMETERS 

3.1 Gas models 

The aims of Sections |3.1| and |3.2| are to test the dependence 
of rh,2/rh on the cluster parameters and compare the results 
with the theory presented in Section [2.1[ The simulations in 
this section ||3.1|) were run using a two-component gas code 
(see |Heggie fc Aarseth| ( 11992J and |Breen fc Heggie|f2012a) ). 
In all cases, the initial conditions used were realisations of 



the Plummer model (Plummer 1911 Heggie fc Hut 20031 



The initial velocity dispersion of both components and the 
initial ratio of density of both components were equal at all 
locations. The choices for the ratio of stellar masses (1712/ m\) 
are 2, 5, 10, 20, 50 and 100. The values of the total mass ratio 
{M2/M1) used in this section are 0.5, 0.1, 0.05, 0.02 and 0.01, 
though the first two values are outside the parameter space 
of interest in most of the present paper. The value of ?'2/rh, 
which was found to be approximately constant during the 
post collapse phase of evolution (see Fig. IT]), was measured 
for the series of models and is given in Table [l] 

The results in Table [l]are plotted in FigTlS] As can be 
seen in Fig. [2] (and Table [l| the results are in qualitative 
agreement with the theory in Section [2. 1[ in the sense that 
for M2/M1 < 0.1 there is an increase in the values of rh,2/ni 
with increasing 7712/7711. For M2/M1 = 0.5 the trend is quali- 
tatively different than for M2/M1 < 0.1; there is an increase 
in the values of rh,2/rh with decreasing 777,2/7711. The theory 




1.5 2 2.5 

Time (t^j, ;) 



Figure 1. Top: r\ (top line) and rj, 2 (bottom line) vs time 
(units trh.i) of g3,s models with A'' = 32fc, 7712/7711 = 10 and 
M2/M1 = 0.02. rji and rij 2 are given in N-body units. Initially 
rjj and rh 2 have the same value, but mass segregation quickly 
decreases rjj 2 , after which it reaches an approximately steady 
value. Bottom: rj^ 2/7'h vs time (units trh,i)- Core collapse (see 
Section|2.1|| occurs at Ri Krh.i; shortly before this rji 2/7"h reaches 
a nearly constant value. 



in Section |2.1| cannot be expected to apply in this regime, 
as the light component does not dominate. In this regime 
the decrease of rh.2/i'h with m2/mi can be explained quali- 
tatively by the fact that if the BH have larger stellar masses 
then there is a stronger tendency towards mass segregation 
(see |Breen fc Heggie| ( [2012a| ) for a discussion of this topic). 

We now consider the comparison with theory more 
quantitatively in the regime M2/M1 < 0.1. In Fig. [2] we 
can see that values of rh,2/fh are also in quantitative agree- 
ment with equation |4] for 7712/1111 > 10 in the sense that 
the power law index is approximately confirmed. However 
''h,2/i'h increases more rapidly then is expected by equation 
|4]for 1112/1111 < 10 and M2/M1 < 0.05. This behaviour is 
possibly explained by equation |6J which predicts a different 
power law index for Spitzer stable systems than in equa- 
tion |4] Indeed for the lowest values of M2/M1 the slope is 
approximately consistent with equation |6] Nevertheless the 
stellar mass ratios of interest are 7712/7111 > 10 as realistic 
ratios for systems containing BH sub-systems would be in 
this range. 

The variation of rh,2/i'h with M2/M1 over a range of dif- 
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Figure 2. The variation of rh,2/^h witli m2/mi. Tlie points rep- 
resent the values of rij 2/''h given in Tablell] The hnes represent 
the expected variation of rh,2/^h with ra^frnx (fitted curv es of 
the form fe(m2/nii)''*) as given by the theory in Section |2.l| The 
Hnes have only been included for M2/M1 < 0.1 as this is where 
the theory is expected to apply. The empirical measured varia- 
tion of ^11,2/^11 with m2/m\ is in good agreement with theory in 
Section [|j] when m2/mi > 10 and M2/M1 < 0.1. See text for 
further details. 



0.25 




0.15 



0.05 



2000 3000 

Time (N-body units) 



5000 



Figure 4. rh,2/^h vs time (in N-body units). N-body runs with 
initial values N = 64fc, M2/M1 = 0.02, m2/mi = 20 (thick Une) 
and m2/m\ = 10 (thin line). Time is set so that core collapse 
occurs at t = for both systems. r\^ 2/''h has been smoothed to 
make the plot clearer. The values of rh,2/''h in the graph are in 
approximate agreement with the results from the two-component 
gas model given in Table fl] 
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Figure 3. The variation of r'h,2/^h with M2/M1. The points 
represent the values of rjj 2/''h for the case 1112/ mi = 10 with N = 
32fc, TV = 64fc and A'^ = 128A:. The line represents the expected 
variation of rij 2/''h with M2/M1 as given by the theory in Section 



o 



ferent A'' is shown in Fig.ls] for the case of m2/mi = 10. For 
the case of A*' = 32fc the variation is less than expected from 
equation B] (see Section [2.1[ ). The variation is approximately 
of the same form, i.e. rh,2/rh oc {K'h/Mi)"' , but a ~ 0.3 for 
the case of 1712/1x11 = 10 and N = 327^, which is less than 
the expected value of a ~ 0.6. (Here we ignore the depen- 
dence on the Coulomb logarithms). The variation of rh,2/rh 
with M2/M1 comes into better agreement with equation [4] 
with increasing A'', with the case of A'^ = 128fc being in good 
agreement with the theory in Section |2.1[ This seems to in- 
dicate that the disagreement is caused by small values of N2 ■ 
One of the assumptions under which equation|4]is derived is 
that re ^ rh,2, but it is possible that re, 2 ~ rh,2 for small 7V2 
as shown in equation [5] In this case one of the assumptions 
underlying the theory is not satisfied. 



3.2 rh,2/rh in N-body runs 



In Section 2.1 it was predicted that rh,2/'"h would increase 



with increasing 7712 /mi in a system with fixed M2/M1 and 
A'^. This prediction will now be tested with direct N-body 
runs (see SectionHl. The initial conditions are realisations of 
the Plummer Model with N = 64fc, M2/M1 = 0.02 and two 
different values of m2/mi (10 and 20). We have compared 
the values of rh,2/rh in Fig. El The value of rh,2/rh is indeed 
larger for m2/mi — 20 than for 1712 /mi = 10 as expected. 
This effect was confirmed using a two-component gas code in 
the previous subsection. However mass is conserved in the 
gas models whereas in the more realistic N-body systems 
mass is lost over time. Therefore we need to ensure that 
we are comparing the values of rh,2/rh for both the runs 
at constant M2/M1. For the two runs in Fig. El mass is lost 
from the BH sub-system at approximately the same rate (see 
Fig 11 bottom), as predicted by the theory in Section 2.5 



E 



The mass loss from the light component is negligible over 
the length of these N-body runs: less than 5% of the light 
component is lost during the entire run. Therefore while 
M2/M1 does decrease with time over the runs in Fig.B] the 
values of M2/M1 for both runs are approximately the same 
at any given time. 

The variation of rh,2/'"h with M2/M1 is shown in Fig. 
5| for the N-body run with A'^ = 64fc, ni2/mi = 10 and 
M2/M1 = 0.02. Initially rh,2/rh = 1 before being quickly 
reduced due to mass segregation. The BH sub-system starts 
producing energy when rh,2/rh reaches ~ 0.1 and M2/M1 
begins to decrease at this point. The variation of rh,2/rh with 
M2/M1 is again less then expected from equation|4| with the 
result indicating a dependence of rh,2/rh oc {M2/Mif'^^ . 
Although this is not in agreement with equation |4] if we ne- 
glect the coulomb logarithm, the results are in good agree- 
ment with the gas model. 
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Figure 5. rii.2/^h vs M2/M1. The solid line represents results 
from an N-body run with initial values A'' = 64A;, M2/M1 = 
0.02 and m^jviix = 10. The results are smoothed to make the 
value of rh,2/''h clearer. The dotted line is the best matching 
curve of the form b( M2/M1)'' (where a and b are constants, the 
best match values being b S! 0.33 and a S! 0.28). The points are 
results from the two-component gas model (see Table ml , where 
M2/M1 is fixed. At the beginning of the run rjj 2/''h = 1 before 
mass segregation rapidly reduces its value, whence the vertical 
line segment in the top right corner. 



Table 2. Variation of \(f>c\ and 02/01 with m^lmx for systems 
with M2/M1 = 0.02. The values are measured in the post-collapse 
phase of evolution when the systems have reached a certain size 
(rh ^ 0.83). 
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3.3 Central potential 



In section pTSl we made the assumption that the main contri- 
bution to the central potential is from the light component. 
In order to test this assumption the central potential and 
the relative contribution of each component to the central 
potential (02/</'i), have been measured in a series of two- 
component gas models with M2/M1 = 0.02. These results 
are presented in Table [2] The values in Table [2] were mea- 
sured once the systems had reached a certain value of Vh 
(rh = 0.83). This was done in order to insure the systems 
had reached a similar point in their evolution (see Fig. |6|. 
The variation in (j)c in Table [2] is only about a factor of 1.2 
for fixed A'^ even though 1712/1x11 varies by a factor of 5. The 
variation in 02/01 is higher, but the values are of the same 
order of magnitude as the estimate in Section |2.5| and in- 
deed in satisfactory agreement, considering that M2/M1 is 
higher here. As the energy gener ation rate per unit mass is 
oc mipi/o"! ( Heggie fc Hut|2003 1, for systems with the same 
value of A'^, M2/M1 and mi, as m2 increases the system can 
produce the required energy at a lower central density, thus 
the central potential is expected to becomes shallower with 
increasing 1x12 1 m\ as seen in Table [2] This is also why two- 
component systems with larger 7712 /rrai are stable against 
gravothermal oscillation to higher values of A*' then for lower 
w.2/mi ( Breen fc Heggie|2012a l. 

These values serve as a rough guide to how the varia- 
tion of 1112 /mi will affect the system. If we consider systems 
with fixed A*', M2 and M2/M1 and use similar reasoning as in 
Section 2.5 we would expect systems with lower m2/m\ to 
last slightly longer than systems with higher 7712 /mi. This is 
because the average energy contribution per binary is depen- 
dent on the depth of the central potential. The deeper the 
central potential the greater the average energy contribution 
per binary will be, and therefore it is expected that for fixed 



:0- 
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Figure 6. Top: |0c| vs rjj for gaseous systems with fixed A'^ 
(128fc) and M2/M1 = 0.02; the different values of m2/m\ are 
10, 20 and 50. Bottom: 02/01 vs r^ for the same systems as in 
the top figure. This plot shows that the contribution of the light 
component to the central potential is dominant. 



M2/M1 a two-component system with ni2/m\ — 20 will lose 
mass slightly faster than a system with m2/mi — 10. 



4 EVOLUTION OF THE BH SUB-SYSTEM: 
DIRECT N-BODY SIMULATIONS 

4.1 Overview 

In order to study BH sub-systems we carried out a number 



of N-body simulations using the NB0DY6 code (Nitadori & 
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Table 3. Parameters used for N-body runs. In all cases the ini- 
tial conditions were realisations of the Plummer Model ( |Plummer| 
|1911||. The values of tj-h,! (the initial half-mass relaxation time) 
are calculated using A = 0.02A'^ for the coulomb logarithm. 



N 



m2/mi M2/M1 tj-h i (N-body units) 



32fc 
64fc 
128A: 



10, 20 

10, 20 

20 



0.01, 0.02 

0.01, 0.02 

0.02 



471 

851 

1552 



Aarseth 20121. Because of the computational cost of large 



A'^ simulations we have mostly limited ourselves to runs with 
N — 32k and 64A:. The total mass ratios used were M2/M1 = 
0.01 and 0.02. The stellar mass ratios used were m2/m\ — 10 
and 20. An additional run with N = 128A;, M2/M1 = 0.02 
and 7712 /jTii = 20 was also carried out to increase the number 
of systems with N2 > W^. These parameters are summarised 
in Table [S] and the results of these runs are given in Table 
[4] The evolution of the fraction of mass remaining in BH 
{M2/M2,i) in all runs is shown in Figs 10 11 and 



12 



First we will discuss the qualitative behaviour of sys- 
tems containing a BH sub-system. We do this by consider- 
ing the case 7712/7711 = 20, M2/M1 = 0.02 and N = 64fc, 
which we refer to henceforth as (20,0.02,64fc). The graphs 
of rh and r^ against time for this run are shown in Fig. W\ 
The BH population quickly segregates to the centre of the 
system causing core collapse to occur in the BH sub-system. 
For the parameters of this model, this takes approximately 
0.3frh,i (where irh,i is the initial half-mass relaxation time) 
to occur. The collapse time for other parameter choices is 
discussed at length in |Breen fc Heggie] ( |2012a[ ). In Fig.lTlthis 
occurs at 270 N-body units. We will refer to this as the first 
collapse. This is followed by a phase of powered expansion. 
This can be seen in Fig. W\ where, after the first collapse, 
both Tc and rh increase up until ~ 5000 N-body units. As 
BH escape, the BH sub-system becomes less efficient at pro- 
ducing energy (see Section |2.3[ ) and the rate of expansion 
decreases. The core stops expanding and begins to contract 
again at ~ 7500 N-body units; at this stage there is only 
15% of the BH sub-system remaining (see Fig. |11[ bottom, 
solid line). Most of the remaining BH escape before « 9500 
N-body units leaving the system with a single remaining BH 
binary from ~ 11500 N-body units. The contraction of the 
core that begins after ~ 7500 N-body units shall be referred 
to as the second core collapse or recoUapse. As with the first 
core collapse the core is contracting because there is not 
enough energy being produced to meet the energy demands 
of the cluster. As the core (which is dominated by the low 
mass stars as most of the BH have escaped) becomes smaller 
towards the end of the run the remaining BH binary starts 
to interact strongly with the light stars, producing energy 
more efficiently. This causes the more rapid increase in r^ 
seen towards the end of the plot. The contraction of the 
core was still ongoing at the end of the run. The core will 
presumably continue to contract until balanced evolution is 
restored. If the last remaining BH binary is providing most 
of the energy, then how long that binary persists in the sys- 
tem depends on the hardness of that binary. Assuming the 
last BH binary is only slightly hard it is possible that the 
core contracts sufficiently for a single BH binary to produce 
the required energy to power the expansion of the system. 
Ultimately the BH binary will become hard enough to cause 
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Figure 7. N-body run (20,0.02,64fc), i.e. with N = 64fc, 
m.2/mi = 20 and M2/M1 = 0.02. Bottom: re vs time (N-body 
units). At t = 270 N-body units core collapse occurs. After most 
of the BH have escaped (t Ri 7500) the core starts to re-collapse. 
Top: rjj vs time (N-body units), rji initially expands rapidly be- 
fore gradually slowing down as the BH sub-system dissolves. From 
t ft! 10000 to ft! 15000 there is little change in rj^ as the system 
is no longer in a balanced energy generating phase of evolution. 
The expansion after t ft; 15000 results from a single remaining 
BH binary which becomes more active as the core collapses. 



its ejection from the system. However if it is extremely hard 
it is likely that the binary gets ejected from the system dur- 
ing the second core collapse. If this happens the collapse of 
the core will continue until light binaries are produced as in 
a one-component model. 



4.2 The rate of loss of BH 

Now we compare the values of M2 with the theory of Sec- 
tion 



2.5 



(equation 10 1. M2 is estimated by calculating the 
average mass loss rate over the time taken (from the start 
of mass loss) for the BH sub-system to lose 50% of its ini- 
tial mass (i.e. 0.5M2,i/Tc,Q%; where Te,Q% is the time taken 
from tec until 50% of the BH sub-system has escaped, see 
Table B] for details). Note that there is a small systematic 
error introduced by measuring mass loss from the point at 
which it first occurs. The values of M2 (in units of lO^^t^^^.^) 
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Table 4. Results from N-body runs with parameters given in Table^ The values given are stellar mass ratio {m2/mi), initial total mass 
ratio (M2/M1), initial total particle number (N), initial number of BH {N2), core collapse time (ice), the time at which 50% (750%) 
and 90% (790%) of the initial BH total mass has escaped, the recoUapse time of the system, the rate of mass loss from the sub-system 
— M2t C (s^6 Section pk and the number of the figure which plots the fraction of remaining BH mass {M2/M2^i) with time. Times for 
TgQojj, Tgo% and the recoUapse are given in N-body units and given in brackets in units of i^ji ;. Times for Tc^o%, Tqq% and the recoUapse 
are measured from the time at which core collapse finishes. The value of Tgo% for the case m2/mi = 20, M2/M1 = 0.01 and N = 32A; 
(marked with *) is actually the point where 88% mass loss occurs; after this point all that remains in the system is a single binary BH. 



The values of M2 are given in units of 10 ^ N-body units, (or 10 ^ t^^ j for the values in brackets); these are measured between the loss 

of the first BH and 50% of the BH, by M2 = —0.5M2^i/Tc^Q%. The values given in the subscript and superscript are the upper and lower 
90% confidence limits assuming that BH escape is a Poisson process. The values of (^ were measured by assuming E/\E\ Si rh/''"h (which 



0.138iV2 d{r,^) , r-^ ,, , , 

(Ri i^ij — N-body units) between 2tcc and ic 



holds if \E\ <x GM^ Jtt^ and M is small) and evaluating 

In A 3 dt 
3 
evaluated by taking the slope of the best fit line to rj^ ; the typical errors with the fitted lines were small, < 3% 



-^50%- —^ was 
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rh,!'' 



are plotted in Fig. [s] The error bars (estimated as stated 
in the caption of Table |4| are large because N2 is relatively 
small, and most data points are con sist ent with a value of 
approximately 3x \Q~^t^^-^. Equation 



10 



with canonical val- 



ues of a, (3 and C, implies that the values in Fig. |8] should b 
nearly 6.1 x 10"" 
the discrepancy. 



and we will now consider reasons for 



First, this estimate can be improved upon by taking 
into account the fact that the system expands as energy 
is being generated, increasing the relaxation time and in 
turn decreasing the mass loss rate. The improved estimate 
can be calculated by using equation [12] to estimate the time 
taken for half the BH to be lost and evaluating M2 as was 
done in Table |4] This results in a slightly smaller estimate 
of5.4x 10~'^i^;, which is still significantly larger than most 
of the values in Table |4] 

Another factor is that the values of C, for most of the 
runs are smaller than the canonical value used for the esti- 
mate (i.e. C, ~ 0.09). Equation 12 predicts an approximately 



linear dependence of M2 on ij, and this is clearly confirmed 
in Fig. [9] The solid line, which represents the predicted val- 
ues of M2 with varying C,, nevertheless lies above all the 
numerical results and outside the confidence intervals for all 
but a few of the runs. However by adjusting the value of 
q//3 (in equations 10 and 12 a and /3 only appear in the 



form a//3) from a//3 « 0.068 (the value estimated on the 
basis of theoretical arguments) to aj fi ~ 0.051 the theory 
comes into very good agreement with the values of M2. This 
can be seen in Fig. [9] where the dashed line represents the 
predicted values of M2 based on a value of a//3 « 0.051. 
The discussion of Section [2.51 makes it clear that the canon- 
ical values of a and, especially, /3 are subject to uncertainty, 
the latter resting entirely on approximate theoretical argu- 



ments. The suggested revision of aj fi cannot be ruled out 
on these grounds. 

Equation [13] allows us to test the theory constructed in 
Section [2.5| in a (^-independent way. In Fig. [l3] the observed 
dependence of M2 on rh is in satisfactory agreement with 
the predictions based on equation |13[ 

The lower values of C, in Table |4] may result from sys- 
tems in which the BH sub-system is incapable of producing 
the required energy for the system to achieve balanced evo- 
lution. This could possibly be due to the small values of A'^2; 
most of these models reach values of N2 (at time T-^q%) be- 
low the point at which the theory of Section [23] is expected 



to apply (see Section 2.31. If a system is not in balanced 
evolution it is expected to undergo contraction of the inner 
Lagrangian radii relative to rh, qualitatively as in conven- 
tional core collapse. This is illustrated in Fig. [M] for three 
of the N-body runs in Table [4] The systems with smaller 
values of C, show greater contraction. Note that the values 
C, are only evaluated over the period to Tc^q% and appear 
to decrease after T<^q% (~ 3000). Indeed this is what would 
be expected as the expansion is affected by the weakening 
energy generation (see also Fig. [7|. 

In this section we have assumed that mass segregation 
concentrates the BH in the centre of the system by the time 
of the first core collapse. Though this prevents further con- 
traction of the central BH sub-system, this is not true of 
all the BH, as discussed in Appendix [B] and IMorsc her et al| 
(|2012|). The outermost BH can continue contracting after 
core collapse has occurred, indicating that mass segregation 
can continue in the outermost parts of a system for a while 
after core collapse. This results in additional heating which 
is not associated with energy production. The effect of this 
heating is expected to be small for the models in Table [3] 
due to the small particle number, although the effect may 
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Figure 8. M2 (units of 10 '^t^.j^ j) versus initial value of 
N2. Error bars indicate confidence limits (see Table W] for 
details). Circles represent {20, M2/M1, N), squares represent 
(10, M2/Mi,Af), filled symbols represent (m2/mi, 0.02, Af), 
unfilled symbols represent (m2/rrti, 0.01, A'^), larger symbols 
(7712/"^!, M2/M1, 64A:) (with the exception of the largest circle 
on the right which corresponds to (20, 0.02, 128fc)) and smaller 
symbols {m2/ra\, M2/M1, 32A;). For cases with the same or simi- 
lar initial values of N2 some of the values were adjusted by < 10% 
to stop the symbols from overlapping. To a first approximation 
M2 is independent of 1712 /mi, M2/M1 and A^. 




Figure 9. M2 (units of 10 '^t^.j^ j) versus value of C,; error bars 
indicate 90% confidence limits (see Table [4] for details). The sym- 
bols represent the same runs as in Fig. Is] the solid line represents 
the predicted values of M2 using the value of a//3 = 0.068 (based 
on theoretical arguments), and the dashed line represents the pre- 
dicted values of M2 using the value of a//3 = 0.051 (the empirical 
value). For cases with the same or similar values of C, some of 
the values of C, were adjusted by ±5% to stop the symbols from 
overlapping. See text for details. 
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Figure 10. Fraction of initial mass remaining {M2/M2 i, where 
A/2,i is the initial mass of the heavy component) vs time (in N- 
body units) for the cases N = 32fc with M2/M1 = 0.01 (Top) 
and M2/M1 = 0.02 (Bottom). In both figures the dashed line 
represents m2/m\ = 10 and the solid line represents m2/m\ = 
20. T = is set as the time when first mass loss occurs, which is 
at approximately the same time as the first core collapse. 



be more significant in larger systems and may be enhanced 
by the presence of a mass spectrum. 

4.3 Lifetime of BH sub-systems 

Now that the dependence of M2 on cluster properties has 
been discussed, it is natural to move on to considering how 
long the BH sub-system lasts. For this purpose we shall de- 
fine the life time of the BH sub-system as the time taken 
from the core collapse of the BH sub-system (which occurs 
at tec) until the BH sub-system has lost 90% of its initial 
mass (Tgo% ) . These values are given in Table B] where it 
can be seen that TgQ% ~ lOirh.i- Equation |12| in Section 
|2.5| can be used to estimate the life time of the sub-system 
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Figure 12. Fraction of initial mass remaining (M^jM^^i, wliere 
M2,i is the initial mass of the heavy component) vs time (in N- 
body units) for the case N = 128fc with M2/M1 = 0.02 and 
m2/mi = 20. T = is set at the time when the first mass loss 
occurs, which is at approximately the same time as the first core 
collapse. 
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Figure 11. Fraction of initial mass remaining {M2/M2 i, where 
M2,i is the initial mass of the heavy component) vs time (in N- 
body units) for the cases TV = 64fc with M2/M1 = 0.01 (Top) 
and M2/M1 = 0.02 (Bottom). In both figures the dashed line 
represents m2/mi = 10 and the solid line represents 7712/ mi = 
20. T = is set at the time when the first mass loss occurs, which 
is at approximately the same time as the first core collapse. 
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(using a/13 = 0.051). For M2/M1 = 0.01 the theory pre- 
dicts Tqo% ~ 2.2ti.h.i and Tgo% ~ 5.2irh,i for M2/M1 = 0.02. 
These values are significantly smaller than the values seen 
in Table HI As stated in the previous subsection most of the 
values of (" in Table|4]are below the value used in Section [2?5l 
Adjusting ( to 0.05 in equation [l2] increases the predicted 
values of Tgo% to 4.0irh,i for M2/M1 = 0.01 and 9.3irh,i for 
M2/M1 = 0.02. These values are still significantly smaller 
than those given in Table 111 with the corresponding value of 
^. The difference between the empirically found values and 
the theoretical estimates might be accounted for by the fact 
that the theory assumes a constant value of ("; however ( is 
expected to decrease as the BH sub-system evaporates (see 
Section 4.11. This behaviour is illustrated by the behaviour 



Figure 13. Evolution of A/2 vs rj, for the N = 32k and 64fc mod- 
els in Table [4I Th e thick dashed line is the theoretical prediction 
(see equation |13| for the initial value of M2/M1 = 0.02 and the 
thick solid line is the theoretical prediction for the initial value 
of M2/M1 = 0.01. The value of rj, ; used in equation |13| was 0.77 
and the empirical value of a//3 S! 0.051 was used for all models. 
In all cases the behaviour of the N-body runs is in approximate 
quantitative agreement with the predicted behaviour until there 
are only a few BH remaining. 



of rh in Fig. (tI There is a hint that the same decrease of 
C, with decreasing A^ may also be present in one-component 
models ( Alexander fc Gieles|2012 l. Also as can be seen from 
the values of 1^ in Table |4| some systems appear to be inca- 
pable of achieving balanced evolution at any time through- 
out the loss of the BH sub-system. 

The expected evolution of a BH sub-system can be il- 



© 2012 RAS, MNRAS 000,[lH20l 



14 P. G. Breen and D. C. Heggie 




Time (N-body units) 

Figure 14. Relative contraction in Lagrangian radii (40%, 30%, 
20% and 10%) over time (N-body units) for (20,0.02,64fc) solid 
line, (10,0.02,64fc) dashed line and (10,0.01,64A:) dot dash line. 
T = in all models is set at the time mass loss starts from the 
BH sub-system. Radii are measured in units of the half-mass ra- 
dius. The values of (J for the three runs are 0.08, 0.05 and 0.03, 
respectively. The parameter (^ measures the dimcnsionloss expan- 
sion rate of the half-mass radius and the figure shovirs that slowr 
expansion is associated with relative contraction of the inner La- 
grangian radii. 



lustrated using Fig. [9] If we assume that the BH sub-system 
is capable of achieving balanced evolution, after the for- 
mation of the BH sub-system ( and —M2 are predicted to 
rapidly reach the balanced evolution values of ^ ~ 0.09 and 
— M2 « 4 X 10~^t~^^^ (just to upper right of the large filled 
circle). As the BH sub-system loses mass it will eventually 
reach the point where it is no longer capable of generating 
the energy needed for balanced evolution. After this the sys- 
tem will move down the dashed line towards the origin. As it 
does so the rate of mass loss decreases, prolonging the life of 
the BH sub-system. This picture may explain the longer life- 
times given in Table [4] and is consistent with the evolution 
of rh in Fig. [7] 

Finally we briefly consider the recoUapse time of these 
systems (see Tablepl. This is the time between the first and 
second core collapse. It can be interpreted as approximately 
the time it takes for the system to achieve balanced evo- 
lution once the BH sub-system has been exhausted. Most 
of the N-body runs do not reach the second core collapse, 
and therefore mostly lower limits on the recoUapse time are 
given in Table |4] From the results in Table HI this time is at 
least roughly the same time as the core collapse time of a 
one component Plummer model (« 15trh,i see [Heggie fc Hut 



2003) but can be longer because of the offsetting effect of 



BH heating. 



5 GRAVOTHERMAL OSCILLATIONS 

The conditions for the onset of gravothermal oscillations in 
two-component models have been studied by Brecn fc Heg-^ 
gie (2012a I, who found that the value of A'^2 (the number of 



heavy stars) could be used as an approximate stability con- 
dition (where the stability boundary is at A'^2 ~ 3000) for a 
wide range of stellar and total mass ratios (2 Sj m2/mi ^ 50 
and 0.1 ^ M2/M1 ^ 1.0). [Breen fc Heggie| ([2012b[), who 



researched the onset of gravothermal oscillation in multi- 
component systems, found that the parameter called the 
effective particle number A^of (defined as M/rrimax) could 
be used as an approximate stability condition for both 
the multi-component systems they studied and the two- 



component models of Breen & Heggie (2012a I. The stability 
boundary they found was at A^cf ~ 10 , which is also con- 
sistent with the stability boundary of the one-component 



model at A = 7000 (Goodman 19871. However both those 



stability conditions relied on the assumption that the heavy 
component (or heavier stars for the multi-component case) 
dominated the evolution of the system, in the sense that the 
heavy component determined the rate of energy generation. 
This is not the case for the systems considered in the present 
paper as the total mass in the heavy component is so small, 
and so A^ef will not be considered further. (In the systems 
considered in this paper the heavy component dominates the 
production of energy, but the light component controls how 
much energy is created.) We will now investigate the onset 
of gravothermal oscillation for the systems of interest in the 
present paper (where M2/M1 <^ 1.0 and m2/mi > 5). 

The critical number of stars (Adit) at which gravother- 
mal oscillations first manifest was found for the gas models 
with M2/M1 = 0.05,0.01 and, m2/mi = 5, 10, 20, and 50 
(see Section |3|. The results are given in Table p] The values 
of A2 at Acrit for the runs in Table [5| are given in Table [6] 



The values for M2/M1 — 0.5 and 0.1 from [Breen fc Heggie 
(2012a I have also been included in Tables [5[ and [6] for ref- 



erence. For fixed m2/mi, the system becomes unstable at 
a roughly fixed value of N2 (iV2 « 1500 for M2/M1 = 0.05 
and A2 ^ 900 for M2/M1 = 0.01). These results suggest 
that N2 still provides an approximate stability condition (for 
fixed M2/M1) even for models in which the heavy compo- 
nent only makes up a tiny fraction of the system. 

Given that the theory in the present paper is built 
around the assumption that the light component determines 
the evolution of the BH sub-system, it may be surpris- 
ing that the appearance of gravothermal oscillations seems 
solely determined (for fixed M2/M1) by the number of stars 
in the heavy component. However this may be explained if 
one considers the unique structure of systems containing a 
BH sub-system. The presence of a BH sub-system tends to 
produce a system with two cores, one small BH core and 
another much larger light core, which is larger than the 
half mass radius of the BH sub-system ([Merritt et al[[2004| 



Mackey et al 2008 



15 in the present paper). 



also see Fig. 

As gravothermal oscillation in a one-component system re- 
quires a small core to half mass ratio ( |Goodman[[1987[ ) , the 
light system itself is expected to be highly stable against 
gravothermal oscillations. However as the BH sub-system 
has to meet the energy generation requirements of the en- 
tire system, the BH sub-system can have a very small ratio 
of core radius to half mass radius (see Section [2|. If the 
onset of gravothermal oscillation is a result of the BH sub- 
system itself becoming unstable then it would be expected 
that the BH sub-systems would have similar structure at 
the stability boundary. This is indeed the case as can be 
seen in Fig. [15[ which shows the post collapse density profile 
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Table 5. Critical values of A'' in units of 10^ 
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Table 6. Values of N2 at A^crit in units of 10^ 
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20 


3.2 


1.8 
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1.6 
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3.6 


2.0 


1.7 


0.8 



of two systems (with m2/ra\ = 50 and ra2/mi — 10) near 
the stabihty boundary (i.e. A*' is slightly smaller than A'^crit): 
the profiles of the heavy component are almost identical (in 
terms of density contrast, i.e. p2/ph,2) whereas the profiles 
of the light component are significantly different. In fact if 
one were to plot the BH sub-systems in units of pc,2 vs rc,2 
the BH sub-systems would be nearly indistinguishable. We 
study this more quantitatively below. 

The Goodman stability parameter ( Goodman] 1993 1 (or 
a somewhat modified version ( Breen fc Heggie|2012a|b l) has 
been found to provide a stability criterion. The Goodman 
stability parameter is defined as 

_ -Etot/irh 

where Ec is the energy of the core. The critical value for the 
one-component model is logjQ e ~ —2. This condition was 
also found to apply for the Spitzer stable two-component 



models studied by Kim, Lee & Goodman (19981. However, 
[Breen fc Heggie| ( |2012a[ ) found the critical value of e to vary 
for the Spitzer unstable models they studied. They found 
that by slightly modifying the definition of e (^2) a much 
improved stability criterion could be found, with a critical 
value logjQ 62 « —1.5. We can test a version of these pa- 
rameters for the BH sub-systems by suitably modifying e 
to 



ebh 



-E'BH/irh,2 
Ec,2/tic,2 



where -Ec,2 and -Bbh are the energy of the BH core and 
the total energy of the BH sub-system respectively, ebh was 
measured for a range of systems with M2/M1 — 0.01 and 
the results are presented in Table IT] along with the values 
of rc,2/»"h,2 (following the discussion of the previous para- 
graph). As can be seen in Table [7] all the systems with large 
enough m2/m\ have similar values of logj^Q ebh and rc,2/'"h,2 
at their corresponding value of A^'crit which supports the as- 
sertion that the onset of gravothermal oscillation depends 
on the structure of the BH sub-system. 

The critical values of log^g cbh are larger than the val- 
ues of log JO e found for two-component models by |Kim, Lee| 
fc Goodman] (|1998|) and [Breen fc Heggie| (|2012a[), and that 



found for the one-component model by Goodman ( 1993 1 by 
approximately 1.7 dex. Also the critical values of rc,2/r-h_,2 
are larger than the corresponding critical value for a single- 
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Figure 15. PostcoUapse density profile in gas models of two- 
component systems with M2/M1 = 0.01 near the onset of 
gravothermal oscillations for m2/mi = 50 {N = 4.3 X 10^, top) 
and 1712/1711 = 10 [N = 8.5 X 10^, bottom). The following is 
shown in the plot: p\ (thin line), p2 (thick line), ptot (dashed 
line), core radius of heavy component rc,2 (x), core radius of 
light component rc,i (*) and rij 2 (-I-). The core radii have been 

defined as r^^i = ^/Ocr^ ^/(47rpc,i). For the case of m2/mi = 50, 
the BH sub-system creates a density hole in the light component; 
the density of lights in the centre is approximately a factor of 2 
less then its highest value (which occurs at log r/rj, ~ —0.5). The 
remarkably large value of r^i for this case results from the low 
value of Pel and a high value of u^ j caused by the presence of 
the BH sub-system. 



component system ( Goodman|1987 1, i.e. Tc/rh ~ 0.02. This 
may be because the maximum radius of the isothermal re- 
gion in the BH sub-system is larger than rh,2, as was hinted 



by Breen fc Heggie ( 2012a I . If the condition for gravothermal 



instability is that the density contrast across the isothermal 
region exceeds some critical value, and if the edge of this re- 
gion is well outside rh,2, then it can be understood why the 
critical value of rc,2/'"h,2 is larger than Goodman's value. To 
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Table 7. Values of logj^Q cbh and ?'c,2/''h,2 near A^crit for systems 
with M2/M1 = 0.01. For the corresponding values of Af^rit see 
Table [5] The results in this table indicate that gravothermal os- 
cillation manifests once a certain value of rc,2/'"h,2 (or logio ^) i^ 
reached. See text for details. 



m2/m\ 
loglO '^BH 



5 
0.14 
-0.21 



10 20 

0.12 0.12 

-0.30 -0.31 



50 
0.11 
-0.32 



investigate this, the size of the isothermal region was mea- 
sured for (50,0.01,4.3 X 10**), which is shown in Fig.[l5|(Top). 
The edge of the isothermal region (riso) was defined as the 
radius at which g\ reaches 80% of it central value. This gives 
rc,2/riso ~ 0.022 which is consistent with the value of rc/r\ 
found by Goodman ( 1987 1 for the one-component model. 
(For a one-component gas model, r^/viso — 0.016 near the 
stability boundary.) 



6 CONCLUSION AND DISCUSSION 

6.1 Summary 

In this paper we have studied systems intended to resem- 
ble those containing a significant population of black holes 
(BH), i.e. two-component systems with one component be- 
ing the BH and the other the rest of the stars in the sys- 
tem. It was argued in Section |2.4| that mass loss by evap- 
oration due to two body relaxation in the BH sub-system 
does not cause significant mass loss of BH and can be ne- 
glected. The principal mechanism for removing BH, for the 
models considered in the present paper, is superelastic en- 
counters involving BH binaries and single BH in the core of 
the BH sub-system. By considering these systems to be in 
balanced evolution, predictions were made regarding the BH 
sub-system, for example the escape rate of BH (see Section 
2.5 1. Some of the potential limitations of the theory were 
also discussed in Section [2.31 

The theory in Sections |2.1| and |2.2| makes predictions 
about the structure of the BH sub-system, particular re- 
garding the variation of rc,2/'"h,2 with m2/m\ and M2/M1 
(see equation^. (Here the subscripts 2 and 1 refer to the BH 
and the other stars, M , m denoted the total and individual 
masses, and re, rh the core and half-mass radii, respectively). 
The theory was tested in Section [3.1 [ with gas models and 
was found to be in good agreement with theory under the 
condition that -012 /mi > 10 and that A'' > 128fc (see Figs [2] 
and^. The disagreement with the theory outside of those 
conditions may be attributable to small A'^2 and the fact that 
the systems become Spitzer stable at low m2/mi. One of the 
assumptions of the theory was that the BH sub-system only 
made a small contribution to the central potential (see Sec- 
tion ! 



2.5 1, which was tested in Section 3.3 by measuring the 



contribution to the central potential of each component in 
a series of simulations. 

It was argued in Section |2.5| that the rate of mass loss 
from the BH sub-system should be approximately indepen- 
dent of the properties of the BH sub-system (i.e. M2/M1 and 
1712/7711). This theory only requires that the light component 
regulates the rate of energy production and does not rely on 
the stronger assumption that energy is transported through 
the BH sub-system as outlined in Section [2. 1[ In Section[4.2| 



the results of a number of N-body runs were presented (see 
Table [s] and Table B| and the results were used to test the 
predicted mass loss rates from Section [2.5| With the excep- 
tion of systems with m2/m\ — 10 and M2/M1 — 0.01, the 
mass loss rates for the BH sub-system were all consistent 
with a value of A'l2 ~ 3.0 x 10~"^i~j^^, where i^h is the half- 
mass relaxation time of the entire system. However most of 
the runs had a lower value of the dimensionless expansion 
rate (^ than expected and M2 was found to vary approxi- 
mately linearly with (" (see Fig. |9|. Once the variation of 
^ was taken into account and the value of q//3 (where a, 
P are dimensionless parameters determining the energy and 



central potential) was adjusted to 0.051 (see Sections 2.5 
and 4.21, there was good agreement between the empirical 
values of M2 in Table |4 and the predicted values of M2 



made using equation |12[ The low values of C, seen in Table 
[4] may result from the small number of BH in these systems 
which possibly results in the inability of the BH sub-system 
to maintain balanced evolution. Larger simulations will be 
required before this explanation can be confirmed. 

In Section [5] we considered gravothermal oscillations in 
systems containing a BH sub-system. This extends the pa- 
rameter space of two-component clusters studied by |Breen| 
|fc Heggie] ( |2012a[ ) to lower values of M2/M1 for m2/mi ^ 5. 
The results in this section imply that the gravothermal insta- 
bility manifests when the BH sub-system reaches a certain 
profile (see Fig. |15| and Table [7|. A version of the Good- 
man stability parameter was also tested for systems with 
M2/M1 = 0.01 and was found to provide an approximate 
stability condition, although the critical value was signifi- 
cantly larger than the value measured for one-component 
models. The difference between the critical values for the BH 
sub-systems and the one-component models may result from 
the fact that the BH sub-system is approximately isother- 
mal to larger radii than rh,2. The ratio between rc,2 and 
the radius of the isothermal region (r[so) was measured for 
a selected model and was found to be consistent with the 
critical value of rc/riao found for a one-component system. 
These results indicate that the onset of gravothermal oscilla- 
tion for systems containing a BH sub-system is determined 
by the properties of the BH sub-system. 



6.2 Astrophyical issues 

In the present paper we have made several simplifying as- 
sumptions. Importantly we have ignored stellar evolution 
and the effect of a mass spectrum. A mass spectrum can in- 
crease the rate of evolution of a system ( Gieles et al||2010 1 , 
which by the theory in Section |2.5| would lead to a faster 
escape rate of BH. On the other hand mass loss via stel- 
lar evolution from the formation of the BH can cause the 
system to expand ( [Mackey et al|[2008[ ) increasing the relax- 
ation time in the system. This in turn would reduce the rate 
of energy generation in the system, which by the theory in 
Section [2. 5| would prolong the life of the BH sub-system. 

Another simplifying assumption was not to consider the 
removal of BH by natal kicks, which if large could signifi- 
cantly reduce the retained BH population. The topic of na- 
tal kicks for black holes is still under debate, so here we will 
only give the topic very general consideration. The ejection 
of BH by natal kicks is itself an energy source, which heats 
the system in qualitatively the same way as a BH ejected by 
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Figure 16. The parameter space of two component systems di- 
vided up into different regions depending on their structure and 
the conditions for the onset of gravothermal oscillations. See text 
and Table Is] for details 



superelectic encounters with binaries. Also if a natal kick is 
not significant enough to remove a BH from the system the 
BH would shed much of the kinetic energy gained from the 
kick to the other stars in the system. This is analogous to 
the results of [Fregeau et al| ( |2009[ ), who found that adding 
natal kicks to white dwarfs was an additional energy source. 

Another topic we have ignored is the presence of more 
than one stellar population in many globular clusters. In the 
typical scenario for the formation of the second generation 
stars in a globular cluster ( [Ventura et al|[2001 1 ejecta from 
asymptotic giant branch stars cools and collects in the centre 
of the cluster. If a BH sub-system is already present at the 
centre this could lead to a significant increase in the stellar 
mass of the BH ( jKrause et al|2012[ [Leigh et al|20l2l ) and by 
the theory in Section [2[ an increase in the total mass of the 
BH sub-system would increase its life time. Even in a single 
population scenario physical collisions can occur between 
BH and other stars in the system ( Giersz et al|[2012 ), and 
this also would increase the total mass in the BH sub-system. 

We have not considered systems which contain an in- 
termediate mass black hole (IMBH) alongside a BH sub- 
system. A recent radio survey by jStrader et al (20121 found 
no evidence of IMBH in the three globular clusters M15, 
M19, and M22. However, there could be other clusters which 
contain both an IMBH and a BH sub-system and this would 
be an interesting topic for future work. 



6.3 Classification of two-component systems 

In Section [5| we considered gravothermal oscillations in sys- 
tems containing a BH sub-system. This extends the param- 
eter space of two-component clusters studied by [Breen fc[ 
Heggie (2012a I to lower values of M2/M1 for 7712 /mi ^ 5. 



But there are rather distinct physical characteristics of the 
systems studied in the two papers, as we have already seen 
in the study of gravothermal oscillations (Section |5|. Here 
we attempt to summarise these ideas. 



Table 8. Summary of the different regions in the parameter space 
of two-component systems. The table states whether a given dy- 
namical process is dominated by the light or heavy component 
of the system, t^ji represents the two-body relaxation time within 
rjj, trc represents the two-body relaxation time within re, GTO 
stands for gravothermal oscillations, with reference to which com- 
ponent contains the large isothermal region which becomes un- 
stable, and the final column represents whether or not the system 
is Spitzer stable; (for Region IV, since N2 ^ 2, Spitzer instability 
is not an appropriate concept). See text for further details. 
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In order to differentiate between two-component sys- 
tems a classification scheme has been devised that divides 
the two-component system parameter space into four regions 
(see Fig. |16[ ). The criteria used to divide the parameter space 
are, (a) whether or not the isothermal region which becomes 
gravothermally unstable is associated with the heavy com- 
ponent (regions II & HI) or the light component (regions I & 
IV), (b) whether or not the system is Spitzer stable (region 
I) or Spitzer unstable (regions II & HI) and (c) whether or 
not the two-body relaxation process within rh is dominated 
by the heavy component (region II) or the light component 
(regions I, HI & IV). The differences between the regions 
are summarised in Table [8| We shall now justify the classifi- 
cation by considering each of the criteria used more closely. 

The simplest distinction to make is between systems 
which are Spitzer stable and these which are Spitzer unsta- 
ble, that is systems which achieve equipartition of kinetic 
energy by mass segregation and those which cannot. [Spitzer[ 
( 1987 1 constructed the following stability condition 



M' 



m2 



,T <0.16f 

Ml \mi 

based on theoretical arguments and some simplifying as- 
sumptions. However a study by [Watters, Joshi fc Rasio[ 
(2000), using Monte Carlo simulations, found Spitzer's con- 
dition to be too strong and suggested a different condition 
of similar form with a different constant and power, for the 
range of stellar mass ratios they studied (m^lrrix < 7). For 
simplicity the Spitzer (19871 condition is used in Fig. 16 
to divide the parameter space. The important differences 
between Spitzer stable and Spitzer unstable systems are, 
first, that equipartition of kinetic energy holds after mass 
segregation (i.e. m2a2 ~ mial) and, second, that for any 
appreciable value of nii/mi the value of M^/M-i has to be 
significantly small. Both the systems in the present paper 
and the systems studied by Breen & Heggie (2012a I fall 



into the more general class of Spitzer unstable systems. Re- 
gion I consists of Spitzer stable two-component systems, and 
occupies the lower left part of Fig. [16[ 

Gravothermal oscillations in Spitzer stable systems were 



studied by Kim, Lee & Goodman ( 1998 1. They argued that 



because of the small values of M2/M1 (and the high values 
of 1712/ mi) the heavy component was confined to the centre 
of the system. They showed that the systems they studied 
became gravothermally unstable once a certain ratio of en- 
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ergy flux at rh and Vc (i.e. e) was reached, and that this 
value was the same as that for a one-component system. As 
the bulk of the system is in the light component this implies 
that the instability results from a large isothermal region 
in the light component. These remarks justify the entries in 
line 1 of Table [SI 

One situation where the concept of Spitzer stability is 
inappropriate is where the number of heavy particles is small 
(A*' ~ 2). As the heavy particles tend to find their way 
to the core of the system and form a binary, the role of 
the heavy component is still significant, because this binary 
becomes the power source for the system lying inside the 
core of the light system. That is why these systems have 
been given their own classification, although they may be 
regarded as the extreme of low M2/M-1 for both Spitzer 
stable and Spitzer unstable systems. Clearly in this case if 
gravothermal oscillations are found, they will result from a 
large isothermal region in the light component, hence the 
entries in line 4 of Table [S] 

Finally the last division in Fig.ITSlis between Regions II 
&i III. The space occupied by Regions II & III consists en- 
tirely of Spitzer unstable models. The difference between 
these two models depends on the value of M2/M1. The 
case where M2 <^ Mi, includes the topic of interest in the 
present paper (i.e. systems containing BH sub-systems). Due 
to small M2/M1 in these systems the light component dom- 
inates at rh and thus the rate of energy generation is regu- 
lated by the light component. In the case where M2 > 0.1 Mi 
the heavy component has a significant effect on the relax- 
ation process within rh, particularly when M2 ~ Mi ( |Breen| 
fc Heggie||2012a L The distinction between the two cases is 
clear when considering extreme values of M2/M1 , but the ex- 
act division between the two is unclear and may have some 
dependence on m2/m\. This is why a shaded area sepa- 
rates the two regions in Fig. |16[ In both cases the onset of 
gravothermal instability is associated with a large isothermal 
region in the heavy component (see Section [s] in the present 
paper and Breen & Heggie (2012al). Another reason for the 
distinction between regions II & III is the theoretical ar- 
guments given in Section |2.6| if we consider tidally limited 
systems, M2/M1 is expected to grow with time for region II 
systems and decrease with time for region III systems. 
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APPENDIX A: ENERGY TRANSPORT IN 
SYSTEMS CONTAINING A BH SUB-SYSTEM 

It has been a fundamental assumption in the present paper 
that, in a system containing a BH sub-system, the major- 
ity of energy generated first flows throughout the BH sub- 
system and is then transferred to the rest of the system via 
two-body relaxation. In this section we will discuss the va- 
lidity of such an assumption. We will do this by considering 
the energy generated by a BH binary as it hardens and is 
ultimately ejected from the system, as was done for the one- 
component case by [Heggie fc Hut ( 2003') . We shall assume 
that the BH binary mostly generates energy by encounters 
with single BH. This is a reasonable assumption as long as 
there are sufficient BH for the central region to be domi- 
nated by them. We can divide the life of the BH binary in 
the system into five energy generating phases as follows: 

(1) After the BH binary is formed the interactions be- 
tween the binary and the other BH will not be energetic 
enough to remove either from the sub-system. During this 
phase all the energy that is generated must be deposited 
within the BH sub-system. 

(2) After phase (1) the binary then starts giving the sin- 
gle BH enough energy to escape the BH sub-system. The 
single BH will typically receive more kinetic energy from an 
encounter with a BH binary than the binary itself. During 
this phase the BH binary remains in the BH sub-system and 
the increase in the kinetic energy of the centre of mass (cm.) 
of the binary is deposited in the BH sub-system. Some of the 
energy of the single BH, however is deposited directly into 
the light component though how much is discussed further 
below. 

(3) At some point the cm. of the BH binary will receive 
enough kinetic energy that it too can escape from the BH 
sub-system along with the single BH. As the density is much 
lower outside the BH sub-system, however, the binary is 
unlikely to deposit much energy until it returns to the higher 
density region in the centre of the BH sub-system. 

(4) The single BH starts receiving enough energy to es- 
cape from the system. The BH binary still escapes from the 
BH sub-system but remains bound to the whole system and 
ultimately returns to the BH sub-system. 

(5) Finally the BH binary escapes from the system and 
the binary contributes no more energy to the system. 



Using the same approach as Heggie & Hut (2003 



page 225, Box 32.1), we will now consider where the energy 
generated by each hard binary is distributed. As stated in 
Section |2.5| the amount of energy generated by each hard 
binary is exp ecte d to be ~ 10m2|(/>c| and we will assume 
that 02 /'/'I 



as in Section 



2.5 



10~ . We will consider in 
turn heating from BH which do not escape the sub-system, 
from subescapers (BH which escape from the sub-system 
but not from the system) and finally from BH which escape 
the system. 



Heating from BH which do not escape the sub-system 
is involved in phases 1 (the single BH and the binary) and 2 
(the binary) . The amount of energy generated during phase 
1 can be estimated in the same way as was done for the 



one component case by Heggie & Hut ( 2003 1 . The result is 



that only about 3% of the total energy generation per hard 
binary is generated during this phase. Including the direct 
heating by the BH binary during phase 2, the total heating 
from BH which do not escape the sub-system is about 5.5%. 

Heating from subescapers is involved in phases 2 (the 
single BH), 3 (the single BH and the binary) and 4 (the bi- 
nary). The total amount of energy contributed by subesca- 
pers is about 49%. Initially a subescaper indirectly heats the 
BH sub-system, by m2<j)2 for a single BH and by 2m2<j>2 in 
the case of the binary. The number of single subescapers is 
expected to be approximately 4.3 and the number of encoun- 
ters which cause the binary to become a subescaper is also 
approximately 4.3. This is because in both cases 4.3 is the 
typical number of encounters needed to increase the binding 
energy of the binary by a factor of 10. This brings the total 
amount of heating (including the direct heating considered 
previously) to the BH sub-system to 18%. 

What happens to a subescaper once it leaves the BH 
sub-system is a point of uncertainty. It will indirectly heat 
the light component up to some maximum radius (which 
depends on its kinetic energy) reached by the BH. It is pos- 
sible that afterwards the subescaper remains on a nearly 
radial orbit and falls back into the BH sub-system, releas- 
ing most of its energy there. In the case of one-component 



systems Spitzer (19871 showed that if a particle is ejected 



from the core its orbit is perturbed by the other stars in 
the system, with the result that the particle misses the core 
at the next pericentre of its orbit. However it is clear that 
the more massive BH will be less significantly perturbed by 
the light stars in the system, and may well return to the 
sub-system at the pericentre of its orbit. If the BH does 
return to the sub-system then most of its energy will be dis- 
tributed there, and the same reasoning will also apply to 
non-escaping BH binaries which are ejected from the core. 
If almost all the heating from subescapers occurs in the BH 
sub-system, then the total heating per binary (including the 
other heating considered previously) to the BH sub-system 
is now « 54%. 

Finally it is fairly easy to estimate the amount of en- 
ergy resulting from escapers, which is 7TI2|(/'c| per single BH 
escaper and 27712 1 0c | for the escape of the binary itself. The 
heating is indirect and heating to each component is propor- 
tional to the contribution of each component to the central 
potential. Per hard binary the percentage of energy gener- 
ated due to escape is the same as for a one-component sys- 
tem, and is approximate 45%: of this about 10% goes into 
the heating of the BH sub-system (as 0c — 0c, 1 — 0.10c). 
This type of heating is involved in phases 4 (the single BH) 
and 5 (the binary and single BH). Therefore the heating 
to the BH sub-system per hard binary may be as much as 
59%. We have assumed (Sections 2.1 and 2.2 1 that all energy 



heats the BH sub-system and is then conducted to the light 
component. Clearly this assumption is only approximately 
valid if the proportion of heating to the BH sub-system (per 
BH binary) is ~ 59%. 
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APPENDIX B: HEATING IN OUTER 
LAGRANGIAN SHELLS 

In this section we shall consider the heating by the BH 
caused by the initial mass segregation. In simulations in 
the present paper the BH are initially spread throughout 
the system with the same velocity distribution as the other 
stars. As the BH are much more massive than the other stars 
the tendency towards equipartition of kinetic energy causes 
the BH to lose kinetic energy to the other stars, which in 
turn causes the BH to fall in the potential well of the clus- 
ter. As stated in Section[2] these systems are usually Spitzer 
unstable, and therefore equipartition of kinetic energy can- 
not be achieved: BH continuously fall in the potential well 
until they are concentrated in the centre of the system. The 
time this process takes depends on the location of the BH, 
as the process depends on the local relaxation time, which 
varies significantly throughout the cluster. For BH which 
start in the outer parts of the system this process takes the 
longest, as the relaxation time is longest there and they have 
to travel the furthest to reach the central region. 

For BH whose orbits lie mostly outside the half-mass 
radius of the whole system, the loss of energy causes an 
increase in the mean kinetic energy, just as for orbits in 
a l/r potential. Also, the relaxation time at the location 
of these BH considerably exceeds the half-mass relaxation 
time. Though the equipartition time scale is smaller by a 
factor of order mi/m2, there is some radius outside the 
half-mass radius at which the equipartion timescale is com- 
parable with the time of core collapse which, for parameters 
of our typical models, can be taken to be roughly O.Strh. 
Throughout core collapse, therefore, we can expect the out- 
ermost Lagrangian radii of the BH component to exhibit 
a steady rise in velocity dispersion, while in intermediate 
Lagrangian shells the velocity dispersion first increases and 
then decreases. In the innermost Lagrangian shells the ve- 
locity dispersion may be expected to decrease throughout 
the time to the first core collapse. 

This behaviour is illustrated in Fig. |B1| with an N- 
body run. This run uses mg/mi = 5, which is smaller than 
would be expected for a BH sub-system; however the small 
value of 1712 /mi allows for a large value of ^^2 (in this case 
N2 = 1000), which more clearly illustrates the behaviour. 
As can be seen in Fig. Bl (Bottom) the squared three di- 
mensional velocity dispersion {V2) initially decreases for the 
Lagrangian radii within ~ r^ and initially increases for the 
larger Lagrangian radii. The heating in the outer Lagrangian 
radii ends roughly when the radii enter rh (see Fig. |Bl| Top), 
after which they show a decrease in v'^ up until t ~ 2500, 
when there is an increase in v^ which is associated with 
core collapse. By the time that core collapse has completed 
(at t ~ 3250) most of the BH (80%) have been segregated 
to within ~ 0.2rh. At this time the outermost Lagrangian 



c 
o 
ft 



radius (90% of Ah) has the highest value of v^ (Fig. [Bl 
Bottom), though part of the increase may be attributable 
to BH ejection. The outermost BH are still undergoing mass 
segregation. Indeed as can be seen in Fig. [bT] (Top) the out- 
ermost Lagrangian radius continuously contracts even after 
core collapse occurs. 

The behaviour discussed in this section is not limited to 
systems containing a BH sub-system; in fact any system with 
a mass spectrum should also exhibit similar behaviour. The 
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Figure Bl. Evolution of the heavy component in a two com- 
ponent N-body run with initial values N = 128A:, m2/mi = 5 
and M2/M1 Ri 0.038. The corresponding mass fractions are 1%, 
2%, 5%, 10%, 20%, 30%, 40%, 50%, 62.5%, 75% and 90% of Ma. 
Bottom: Mean square velocity dispersion inside the Lagrangian 
shell. The 90% shell becomes the hottest shell at t ~ 1450 and 
remains hotter than the next inner shell until well after tec- Top: 
plot of the radii of the Lagrangian shells in the heavy component. 
The dotted line is rj, , the half mass radius of the whole system. 



important point is that continued mass segregation serves 
to heat the low-mass component, and this kind of heating 
may lower the required rate of energy generation from the 
core, as energy is directly injected into the region near Vh- 
This effect has not been taken into account in the theory in 
Section [2] but as the values of N2 are so small in Section [2] 
it is unlikely that this effect is significant. 

This paper has been typeset from a TfTJX/ M^rjX file prepared 
by the author. 
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